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We present a simple calculational scheme for superconducting properties under magnetic fields. A 
combination of an approximate analytic solution with a free energy functional in the quasiclassical 
theory provides a wide use formalism for spatial-averaged thermodynamic properties, and requires 
a little numerical computation. The theory covers multiband superconductors with various set of 
singlet and unitary triplet pairings in the presence of an impurity scattering. It is also applicable to 
analyze experimental results in a rotating magnetic field with help of band structure calculations. 
We demonstrate the application to s-wave, d I .2_ !/ 2-wave and two-band s-wave pairings, and discuss 
the validity of the theory comparing with previous numerical studies. 
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I. INTRODUCTION 

A magnetic field has been a good probe to investi- 
gate a gap symmetry of superconductivity. In partic- 
ular, the low-energy excitations in the mixed state ex- 
hibit unconventional behaviors inherent in different na- 
ture of a vortex core, nodal structure and quasiparticle 
transfer between vorticesiSi&i. In a recent development 
in multiband superconductors'^* 6 * 7 *^, magnetic fields 
play an important role, e.g., weak magnetic fields sup- 
press strongly smaller gap of the multiband supercon- 
ductivity leading to drastic enhancement in field depen- 
dence of the residual density of states*fi*ii, the thermal 
conductivity^ 2 ** 3 - and so on. Moreover oriented field mea- 
surements have performed in order to determine directly 
superconducting gap structuro 14 i 15 i 16 i 17 i 18 i 19 i 20 i 21 i 22 i 23 . 

From the theoretical side various approaches have 
been made to investigate bulk properties under mag- 
netic fields. The Ginzburg-Landau theory is powerful 
but is restricted to regions near the critical fields. The 
Volovik theory is used frequently, in which the quasi- 
particle excitation energy is approximated by taking 
into account the Doppler shift of the local supercon- 
ducting currenl^* 2 **!** 2 ** 2 *^. This method neglects the 
scattering by vortex cores and the overlap of the core 
states, which makes the theory to be valid only at low- 
temperature and low-field regions*^. Meanwhile, fully 
numerical calculations for the quasiclassical transport- 
like equation or the Bogoliubov-de Gennes equation have 
been performed* 2 ** 3 ** 3 ** 3 ** 3 *^. Although these approaches 
are necessary to discuss local structure of the core states, 
it is not easily applicable for transport properties in the 
linear-response theory, impurity effect and extension to 
multiband superconductors with detailed band structure 
because of intractable aspect of numerical computation. 

Motivated by these situation, it is useful to arrange a 
semi-analytic approach to analyze various experimental 
results and to check their consistency on an equal foot- 
ing. For weak-coupling superconductors, there exists a 
reliable approximate analytic solutionis in the quasi- 
classical formalism 36,37 with a certain kind of mean-field 



approach. Although required conditions seem to be valid 
near the upper critical field H C 2, it is meaningful to ex- 
trapolate the theory to lower fields. With understanding 
of its validity, it becomes very useful to discuss spatial- 
averaged thermodynamic properties of superconductiv- 
ity under magnetic fields. The theory covers multiband 
superconductors with various set of singlet and unitary 
triplet pairings in the presence of an impurity scattering. 
The realistic band structure can be included in arbitrary 
direction of magnetic fields. Moreover it is also applica- 
ble to investigate linear-response coefficients in a similar 
manner* 3 ** 3 *********^, which will be presented in a subse- 
quent paper. 

The paper is organized as follows. In the next section, 
we introduce the quasiclassical equations in a general 
form, then we give explicit expressions of the transport- 
like equation and the free energy functional for singlet 
and unitary triplet pairings. In §3 we explain the approx- 
imate analytic solution for quasiclassical equations in the 
presence of the impurity scattering. Then we express use- 
ful formulas for basic thermodynamic quantities. In §4 
we extend the theory to multiband superconductors and 
present a set of formulas in terms of dimensionless param- 
eters for the convenience of the users. The application 
of the method to s-wave, d x 2_ y 2 and two-band s-wave 
pairings and the comparison with previous numerical so- 
lutions are given in §5. The last section summarizes the 
paper. 

II. QUASICLASSICAL EQUATIONS 

The microscopic Green's function contains all the in- 
formation about the single-particle properties. In most 
cases, however, states far from the Fermi energy or 
rapid oscillations inherent in the Fermi surface play lit- 
tle role to superconducting properties. The key idea of 
the quasiclassical theory is that those irrelevant states 
are integrated out from the beginning in the Gor'kov 
formalism 3 ** 3 ** 3 ****^. Then, the basic equations of the 
quasiclassical theory give a complete description of the 
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superconducting properties on a coarse-grained scale. In 
this section, we first review the derivation of the standard 
quasiclassical equations under magnetic fields, and we in- 
troduce necessary notations through this step. Next, we 
restrict ourselves to the cases of the singlet and the uni- 
tary triplet pairings, and then we construct a free energy 
functional as a generating functional for the quasiclassi- 
cal formalism. 



A. Derivation of quasiclassical equations 

Let us start with the mean-field (MF) Hamiltonian in 
the 4x4 Nambu-Gor'kov space with a spatial inhomo- 
geneity (c = h = k# = 1 hereafter) , 

H u ?=\ J dx 1 dx 2 ^(x 1 )^k 12 + A 12 ^jp z ^(x 2 ), (1) 

where the creation-annihilation field operators 
at the position x s are composed as \E^(a; s ) = 
[ip\(x s ),ip\(x s ),ip 1 (x s ),ilj l (x s )}. An operator with 
a hat acts on the Nambu-Gor'kov space, and p (<r) is 
the vector composed of Pauli matrices on the 2x2 
particle-hole (spin) space remaining the other degrees of 
freedom unchanged. The unit matrix is denoted by 1. 
The kinetic-energy matrix is then given by 



K 12 = 5( Xl - x 2 )£ {-iV 2 p z - eA 2 l) , 



(2) 



with a single-particle energy operator £(p), where A s 
is the vector potential at x s and e < is the electron 
charge. The superconducting gap matrix is given by 



A 12 = 



A(xi,x 3 ) 

-tf(x2,Xi) 



(3) 



and its component is defined in terms of the two-body 
interaction V a /3 a >pr(xi, x 2 ) as 

A a p(x 1} x 2 ) = - V al3a >f3>{xi,x 2 ){ij}p l (x 2 )ip al (x 1 )). 

a'P' 

(4) 

Corresponding to the MF Hamiltonian, the Nambu- 
Gor'kov equation reads 



dy 



-^(n - r)p z ] p z G(y, x 2 ; Tl - r 2 ) = 15{x x - x 2 ), 

(5) 

for the imaginary-time Green's function 

G{x u x 2 ;n-T2) = -{T T V(x 1 ,n)&(x2,T 2 )), (6) 

where we have used x s = (x s ,T a ) and y = (y, r). Here 
the operator T T arranges the field operators in ascend- 
ing order of the imaginary time < r < 1/T. The 



self-energy £ has been taken into account to treat the 
impurity scattering in subsequent discussion. 

In order to integrate out the rapid oscillations, we in- 
troduce the center-of-mass coordinate, R = (x\ + x 2 )/2 
and the relative coordinate, r = X\ — x 2l and perform 
the Fourier transformation in the latter according to 



f R (k:iuj n ) = J dr J drf {x ll x 2 



t) e 



-i(k-r—uj n r) 



where Lu n — (2n + l)irT is the fermionic Matsubara 
frequency. A product of functions, b} 2 (ji — r 2 ) = 
J dy J drf ly (Ti — T)g v2 (r — t 2 ), is then transformed into 
so-called "circle-product"—, 



h R (k;iuj n ) = exp 



(V fe > ■ Vr- V fe • VhO 



x f R (k\iu) n )g R >(k'\iuj n ) 



(8) 



k'=k.R'=R 



In translationally invariant systems, the internal mag- 
netic field, B = V x A is only the source of slowly 
varying field dependence. In the lowest order in V.r, 
which is necessary in this paper, the circle-product is 
given simply by the product of each Fourier transforma- 
tions, h R (k) ~ f R (k)g R (k). Noting the transformation, 

£(-iV s p z - eA s ) -> 

£fe + l -v(k) ■ l(-l) s V R + 2ieA R p z \ , (9) 

in the leading order, we obtain the Nambu-Gor'kov equa- 
tion as 



iu n + ev(k) ■ A R j p z - - ^v(k) ■ V R j 1 

A R (fc) - cr R {k- ioj n )] p z G R (k; iuj n ) = 1, (10) 



where k = k-p/\k-p\ is the unit wave vector at the Fermi 
surface and v(k) = Vfc£(fcp) is the Fermi velocity. Here 
we have considered that the self-energy in question has 
weak momentum dependence with respect to fcp, hence 
we have denoted as a R (k;iuj n ) — E(fc; iuj n )\k = k F p z . 

To obtain a closed quasiclassical theory, we must con- 
struct equations for a quasiclassical propagator, g ~ 
J d£G. However, a simple integration of eq. ljTU|l results 
in unphysical divergences due to ^G term in the left- 
hand side and 1 term in the right-hand side. We get 
rid of this inconvenience with the help of the "left-right 
subtraction trick" , which transforms the Dyson's equa- 
tion into transport-like equations for g. The trick is to 
eliminate those divergent terms by subtracting the right- 
hand Dyson equation, in which operators in the brace of 
cq. (JSJl act on the second argument of G from the right. A 
similar argument leading to eq. (|10|l gives the right-hand 
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equation, 

p z G R (k;iuj n ) (i 



& + -v(k) • V fl ) 1 - Ajj(fc) - <7fi(fc; iu n ) 



ev(k) ■ A R ) p z 



1. 

(11) 



The unphysical terms then cancel, and the terms left can 
be ^-integrated to yield the following transport-like equa- 
tions, 

(iu) n + ev(k) ■ A R j p z - A R (k) - a R {k; iu n ) , 

g R (k;iuj n ) + iv(k) ■ V R g R (k; iu>„) = 0, (12) 
where we have defined the quasiclassical propagators as 



g R (k;iuj n ) 



'-p z G R {k\iuj n ). 



(13) 



At the subtraction step, a normalization of g is lost. The 
appropriate normalization is known as g R (k;iuj n ) = —1, 
which can be confirmed explicitly in a homogeneous case 
without the internal field Br£&. 

To complete the quasiclassical formalism, we give self- 
consistent equations for the gap and the superconducting 
current as, 

Aj^(fc) =-tcN tJ2Y1 

n a' (3' 

x (v aPa/0/ (k,k')g^ all3 ,(k';iuj n ))^ , (14) 



Vfl X {B R - H R ) = 

4en 2 N TY,J2( v ^9R aa (k^u; n )) , (15) 



where Hr is the external magnetic field and g % i is un- 
derstood as the (i, j)-element in the particle-hole space. 
The bracket (•••) = J dkv^ 1 (k) ■ ■ ■ / J dkv~ l {k) repre- 
sents the angular average over the Fermi surface and N 
is the density of states (DOS) per spin at the Fermi en- 
ergy. 



B. Cases of singlet and unitary triplet pairings and 
generating functional of quasiclassical theory 

In this subsection, we give only expressions for a uni- 
tary triplet pairing. For a singlet pairing, vector quanti- 
ties should be replaced by scalar ones and a unit matrix 
is used instead of er. 

In the case of a unitary triplet pairing, the spin struc- 
ture of the gap is decomposed as 



Afl Q(3 (fc) = Afl(fc) ■ {(Ti(T V ) a p, 



(16) 



with i(A x A*) : 
also has the form, 



Oaa. The quasiclassical propagator 



9R 



-igR.Sa/3 /h • (<ricr y )ap 
f R ■ (o-io-fjL ig R S a /3 



(17) 



By definitions © and (|13fl . the normal and the anoma- 
lous propagators have the symmetry^, 

g R {k; iw n ) = -g R {-k; -w n ) = 9* R {-k; wj n ), (18) 
and 

f R (k;w n ) = TfR.(-k-iu n ) = Tf f R{-k;iuj n ), (19) 

where the upper (lower) sign is applied for the triplet 
(singlet) state. Note that the gap function satisfies 
Afl(fc) = ^Afl(— k) similar to (|T§|> . According to these 
relations, we consider only oj n > hereafter without loss 
of generality. 

At this point we explicitly take into account an effect 
of impurity scattering. Since an impurity potential is 
short range as much shorter than length scale of R de- 
pendence, we merely follow the standard i-matrix theory 
to treat non-magnetic impurities^^Siii. Decomposing 
the impurity self-energy similar to eq. (|17|) . 



we obtain 



T R 



(n), 
T R M 



<r^ ) • (cria v ) a p 



T R 



{9R) 



2t o cos 2 S + Dr sin 5 



(T R '{lU) n ) = 



R) 



2t cos 2 S + D R sin 6 



(20) 

(21) 
(22) 



with D R = (g R ) 2 + (f R ) ■ (f R ), where r is the quasipar- 
ticle lifetime in the normal state and 5 is the impurity 
phase shift (e.g. 5 = corresponds to the Born approx- 
imation and S = 7r/2 is the unitarity limit). Here a R 

and <t^ r satisfy the symmetry relation similar to i|18|) 
and 119(1 . respectively. For convenience, we introduce 
the renormalized frequency and gap function as 



A R (k;iu n ) = A R (k) -+ 



T R ' 



(23) 
(24) 



then the self-energies can be absorbed formally into the 
"tilde" quantities in quasiclassical equations. 
The explicit expression of eq. (|12|) is then given by 

= 5rA r (A;; iuj n ), C-f R = g R A R (k; iuj n ), 
v(k) ■ V R g R = A R (k; iuj n ) ■ f R - f R - A R (k; iu n ), 

(25) 

with the normalization, g R + fR ■ f R = 1, where we have 
introduced 



C±(k,R; iu> n ) 



± -v{k) ■ [Vfl T 2ieA R ] , (26) 
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and £± — C±(u) n — > ^ur)- Note that one of three equa- 
tions in (|25[1 is independent due to the symmetry rela- 
tions ljl8[) and Ijl9|) and the normalization condition. In 
the case of Br = 0, the solutions of eq. H25[) are given by 



and 



A(fc;iw n ).At(fc; 



1 -1/2 



(27) 



f= A{k : lUJn) g, /t= At ( fc '^") g , (28) 



The interaction leading to the triplet (singlet) pairing 
has the form 



W(M') = 4/(M'), 



(29) 



where P a/3a >pi = (Saa'$p/3 l ±S a /3'5i3 a i)/2 is the projection 
operator to the triplet (singlet) state. The interaction for 
the triplet (singlet) pairing is odd (even) against k — > k 
or k' — » —A;'. Then we classify the interaction and the 
gap function by irreducible representations of the point- 
group^. Considering a certain irreducible representation 
r with the highest transition temperature, T c , we can 
write V(k, k) in a separable form, 



V(k,k') = -VT l ^(kWX(k% 



7 



(30) 



with V > 0, where the basis functions satisfy the or- 
thonormal relation, (^ip* (A;)<y9 7 < (k)^ = <5 77 /. Suppose 

Ait(fc) = A.Rip~f(k), the gap equation is rewritten as 



(31) 



V-'Ar = 7tTN ^ Mk)f R (k;iu>. 



Equations i|15|) , l|25|l and l|31|) constitute the quasiclas- 
sical formalism to describe superconducting states under 
magnetic fields. Alternatively the system of equations 
can be derived as the saddle point of the following gen- 
crating functional per unit volume (measured from the 
energy in the normal state)^, 



A R ,A R ,fR,f R ,A R = 

fdR \ {Br - Hr)2 +V-'\Ar\ 2 

-2irTN ((infc *w„)) + l' R (iu„) 

n>0 ^ 



(32) 



where 



I R (k;iLu n ) = <p;(k)A* R fR + f R - ARcp^k) 



1 



1 + ffR 



fi-C+fR + fn-C-fl] (33) 



1 



2r sin^ 6 



In [cos 2 6 + D R sin 2 6]. 



(34) 



In eq. (|32|l the summation over the Matsubara frequen- 
cies must be cut off at very large but a finite frequency. 
This cut-off procedure can be avoided by a prescription, 



V~ 



No 



ln|£ 



2tfV — 



(35) 



namely, the cut-off is absorbed in T c and the second term 
cancels the contribution from the third term in eq. I|32|) 
at n — > oo. Note that T c is the transition temperature at 
H = without the impurity scattering. 

If the quasiclassical equations, eq. are solved sep- 
arately for given Ar, A r and Ar, eq. (|3*3^l becomes 



7 R = A* R (k) fR + f R - A R (k) 



2^(1 -gn) 



9r 



9R 



A R (k;iw n )-f R + f R - A R (k;iu n ) , (36) 



and fisN coincides with the physical free-energy func- 
tional of Ar, A* r and Ar. In the next section, we show 
that a simple approximation provides an approximate an- 
alytic solution of eq. (|25|) . 



III. APPROXIMATE ANALYTIC SOLUTION 

In the previous section, we have derived the set of the 
quasiclassical equations as the saddle point of the free 
energy functional. Basically, a minimum of the free en- 
ergy functional provides a complete description for ther- 
modynamic properties of superconducting states under 
magnetic fields. In this section, however, we show that 
a simple but a reliable approximation drastically simpli- 
fies the problem. In the following subsection, we show 
that a kind of "mean-field" approach leads to an approx- 
imate analytic solution for quasiclassical equations 12511 . 
Then we discuss self-consistent equations for the impu- 
rity scattering self-energy. Lastly, expressions of various 
thermodynamic quantities are given. 



A. BPT approximation 

The approximation used here was first proposed in 
the Gor'kov formalism by Brandt, Pesch and Tewordt 
(BPT) 34 , and later Pesch obtained analytic solutions in 
the quasiclassical formalism 35 . Although BPT initially 
aimed to describe superconducting states near the up- 
per critical field H C 2, it is meaningful to extrapolate the 
method to lower fields. In this direction, Dahm and co- 
workers have made a comparison with the numerical so- 
lution of quasiclassical equations 2 ^. However, in numeri- 
cal calculation they assumed the structure of the vortex 
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cores and did not solve whole equations self-consistently, 
hence a reasonable comparison of the BPT approxima- 
tion with full numerical solutions is still missing. After 
the approximation and the formulation are introduced, 
we will discuss this point. 

In the BPT theory the internal field Br is approxi- 
mated by its spatial average B. Hereafter a spatial av- 
erage of a quantity Xr is denoted by X or Xr. An 
Abrikosov solutionis is then used for vortex lattice struc- 
tures for all field range above the lower critical field. 
When the internal magnetic field B sets parallel to z 
axis, we write the Abrikosov solution with the Landau 
gauge Ar = (0, Bx, 0) as 



-«'$o[(/-AV)/A], (37) 



where p n = 2im/j3, (3 being the lattice constant in the 
y direction and A = {2\e\B)^ 1 / 2 is of the order of the 
lattice spacing of the vortex lattice. We have taken into 
account anisotropy of the Fermi velocity in the x-y plane 
by transferring the coordinates as 



k 



y 



?eff 



(38) 



where (i = x, y) is the z-component of the coherence 
length proportional to the square-root of (vf(k)) and 
£eff = \/S,x£,y sets equal to A. Here $oN = exp(— x 2 /2) 
is the lowest eigenfunction of a harmonic oscillator and 
the periodicity of the coefficients C„ specifies the type of 
vortex lattice. Since k dependence of the gap function 
is irrelevant in the following discussion, we omit ip 1 {k). 
For the vector gap, the argument below can be applied 
separately to each component. 

In addition to the above approximations, we expect 
that the uniform component g — y~R could describe a 
global suppression of the gap amplitude as an average. 
This is because the higher Fourier K components of gR 
decrease rapidly as exp(— A 2 AT 2 ) and the normal propa- 
gator gR does not contain an important phase variation 
due to vortices^i. On the other hand, the exact spatial 
form of Ar due to its phase winding around vortices 
will be taken into account in determining the anomalous 
propagator f^. 

Within the above approximation, we can solve ana- 
lytically the quasiclassical equations (|25|) . It is suffi- 
cient to solve the first two equations in The argu- 
ment here essentially follows the discussion by Houghton 
and Vekhter for s-wave superconductors with the oper- 
ator formalism^!. Within the above approximations the 
anomalous propagator is given by 



f R = g^(k)7 ] ( a >(iu Jn )£- 1 AR 



(39) 



Here we have assumed that the renormalized gap function 
has the same R dependence as the bare gap Ar and the 
renormalization factors of the frequency and the gap are 
independent of R, i.e., 



(40) 



To calculate the inverse of C+ , we use an operator iden- 
tity for uj n > 0, 



/>oo 

C+ 1 = / dtcxp[-C+i\. 
J o 



(41) 



Then we need to know how the exponential operator acts 
on the gap. 

Since the Abrikosov lattice is a superposition of $o 
at different vortex cores, we introduce the lowering and 
raising operators, 



A f ^ 



A 

71 



(42) 
(43) 



where 7r^, = — iVj' — 2eAi> — ^7Ti/A is the dynamical 
(gauge-invariant) momentum for the Cooper pair. The 
operators satisfy the usual bosonic commutation relation 
[a, a'} — 1. Since ciAr — 0, the Abrikosov lattice l|37|) 
can be regarded as the vacuum state |0) in the language of 
the operator formalism. Applying (aJ) m to the vacuum 
state we obtain the excited states, 



C„e" 



'$ m [(x'-A 2 Pn )/A], (44) 



which corresponds to the m-th excited oscillator states 
centered on each vortices. Similarly, we introduce the 
raising and lowering operators as — (a^)* and b = (a)* 
acting on the conjugate vacuum state, = (0|, since 
(a)*A* R = holds. Owing to the factor e tp " y ' , a spa- 
tial average ensures an orthogonal relation for different 
excited states as 



(m\m>) = |A| 2 (5„ 



(45) 



since only functions <& m and centered on the same 
vortex contribute to the integral. Here we have defined 
the averaged gap magnitude as |A| 2 = |Ah| 2 - Without 
loss of generality, we can take A real, because A appears 
only in the form |A| 2 . Note that the orthogonality does 
not hold without the spatial average. 
In the operator formalism denoting 



v(k) = v±(k) COS0U, v±(k) sin^, vn (k) 
we write 



V2 \ 2.\ 

where we have defined 



(46) 



(47) 



v±(k) — v±(k)y x 1 / 2 cos 2 (p v + x 1 / 2 sin 4> v 



tan< 



X 1/2 tan cp v = x 1/2v y/ v x, 



(48) 
(49) 
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with x = ( v x)/( v y)- Note that v±(k) is of the order of 
[(vDiVy)} 1 / 4 , and we obtain v±(k) = V±(k) and </>' = <j) v 
in the absence of the anisotropy, i.e., x — 1- Using Taylor 
expansions in a* and a, we obtain 



c^Ar = 



2A 



E 



-im<p 



. v/2 



W {m \iu n )\m), (50) 



with u„(fe;iw n ) = [2A/u_i_(fe)]u) n > where W(z) = 
e _z erfc(— iz) is the Faddeeva function and W^ m \z) de- 
notes its m-th derivative. The conjugate version CZ 1 A R 
is obtained by the replacements and \m) — ► (m| 

in eq. (|50(l . Then we obtain 

(51) 



where we have used the formula, 

oo 

— ( ; 



^ m! V 5 



m=0 



W (m) (m„) 



=W«(™„), (52) 



and we have abbreviated A</? 7 (fc) to A(fe). 

Using the above result together with the normaliza- 
tion, g 2 + f R ■ fn = 1 , we finally obtain the normal 
propagator, 



<-/ 



1 



2;V 



j/°>y a )|A(fc)| a W (1) (ifin) 



The quantity I = Ir is also obtained as 



(53) 



7 = 



2A 



|A(fc)| 2 VF(w„) 5 



1 + 9 



+ 2 Wn (»jW-l)(l- ff ). (54) 



In the clean limit v = 1/2toT c < 1, we have 
2g ( 2A 
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1 



Vj_(fc) 



A(fc)| 2 IT(iu„). (55) 



In the absence of B using W{z) ~ i/^z and W {1 \z) ~ 
—ij^pKZ 1 for |z| ^> 1, we recover the uniform solution of 
eqs. <|27[) and l|28|l . On the other hand, we get g = 1 and 
7 = 7' = in the normal state, A = 0. 



B. Self-consistent equations for renormalization 
factors, rf n ' (iuj n ) and r}^ a \iuj n ) 

To complete the system of equations, we need self- 
consistencies for rj^ n \iu> n ) and rj^(iuj n ). In the case 



of non-s-wave pairings, we neglect the renormalization 
of the gap, i.e. (Jr) ~ 0, because the angular average 
with the factor <p 7 {k) gives small contribution to the self- 
energy (the angular average completely vanishes in the 
absence of magnetic fields) . We then obtain 



t/"> = 1 + 



(<?> 



2uj n T cos 2 S + (g) 2 sin S 



7 W = 1. (56) 



On the other hand, in the case of the isotropic s-wave 
singlet <p~/{k) — 1, we consider only the Born limit 5 = 0. 
Retaining the term proportional to Ar (m = compo- 
nent) in eq. I|5U|) . namely, 



2A 



W(iu n )A R , (57) 



(this is exact for an isotropic Fermi surface), we get 



ZUl n T 



(58) 



Note that the Anderson theorem (jy(™)/^( a ) = 1) is re- 
stored in the case of 77 = 0. Thus in all cases in the 
limit 77 = 0, we recover all the results of the standard 
theory for dirty superconductors. In the clean limit we 
have rj( n \iu! n ) = rj( a \ioj n ) = 1. 



C. Thermodynamic quantities 

With the approximate analytic solutions obtained in 
the previous subsection, the free energy ll.'i2[| can be re- 
garded as a function of A (can be taken as real) and 77 for 
given temperature T and the external field H. To obtain 
an equilibrium state, we minimize the free energy i]-i2l) 
with respect to (77, A). In searching the minimum point, 
we need to compute the free energy for given (77, A). To 
do this we first solve the self-consistent equations for the 
impurity renormalization rj^ n \iuj n ) and r]^(iuj n ) at pos- 
itive u>nS, i.e. eqs. ()53f) and l|5ti|) for non-s-wave pairings, 
or eqs. <|53|) and l|58|) for s-wave singlet. This can eas- 
ily be done by an iteration starting from the clean limit, 
= jjia) _ ^ Once we determine rf n > and rf a \ we can 
compute the free energy <|32l) via eqs. (|53|> and 1)540. Note 
that all the numerical calculations can be done with real 
quantities in the Matsubara formalism, since 
is real for real x. 

Once we obtain the equilibrium values of B and A, we 
can compute various thermodynamic quantities as the 
spatial average. The magnetization is given by 



-4ttA7(T, H) = H- B(T, H) 



(59) 
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Viewing the one-body nature of the quasiparticles in the 
present formalism, we can write down the entropy, 

S(T,H) = -4 j°° da;N(cj;T,H) /(^)ln/(|) 



(60) 



with f(x) — (l + e*) 1 via the DOS in the superconduct- 
ing state, 



N(u;;T,H) 
N 



= Re ( g(k; iu n — ► cj + iS) ) , 



(61) 



where we have done the analytic continuation from iuj n 
in the upper half plane to lu + iS, S being an infinitesimal 
positive number. Then the specific heat is calculated 
numerically by C(T, H)/T = dS/dT. 

It is useful to discuss the region B ~ H ~ H c2 in the 
clean limit, where A is suppressed. Expanding the free 
energy in A, we obtain the free energy under magnetic 
fields, 



SN 



(B-H) 



8tt 



-N a(T,B)\A\ 2 + ^(3(T,B)(\Af) 2 . 



2T?' 



where 
a(T,B) 



(62) 



v 7 n— v 



and 

oo 

{3(T,B) = -in 2 TT?J2( 



W{iu n ) 



2A 



(63) 



„=o \ \ v±(k)J 
x W(iu n )W^(iu n )\ip 7 (k)\ 4 



(64) 



The temperature dependence of the upper critical field is 
determined by 



a(T,H c2 ) = 0. 



(65) 



In the limit H c2 — > we recover the BCS results, i.e. 
a = ln(T/T c ) and /? = [7C(3)/87r 2 ](|<^ 7 (fc)| 4 ), where ((n) 
is the Riemann zeta function. The field dependence of 
the magnetization and the gap near H c2 are given as 



AttM = 47ra4 2 (T)7V |A| 2 = 



^a%{T)N T? 



-1 (H c2 -H), (66) 



where a' c2 (T) = da{T, B) / dB\ B =H c2 and (3 c2 (T) = 
P(T,H c2 ). At T — we can perform the summation 
over uj n in i|63[l and (|64|l . Then we obtain 



(67) 



with the upper critical field at T = 0, 



5±(fe) 



where w without (k) denotes the angular average of |u(fe)| 
and 7 = 0.577216 is the Euler's constant, and 



0(0, B) 




|^ 7 (k)| 4 



2\e\B V v ) \ \v ± (k) 
In the limit \ -C 1 using v±/vx ~ X 1 ^ 4 , we obtain 
ff c2 (T = 0; X ) 



tf c2 (T = 0;l) 



-1/2 _ 




(69) 



(70) 



IV. EXTENSION TO MULTIBAND 
SUPERCONDUCTORS 



In a system having plural Fermi surfaces, each band 
could have different symmetry and/or magnitude of gaps 
due to symmetry reason^2^Si^i. Since those Fermi sur- 
faces have different shape in general, the Cooper pair 
may be formed among the same band in order to reduce 
the energy associated with the center-of-mass motion. In 
this case the MF Hamiltonian QJ is diagonal in the band 
indices £, and the quasiclassical propagators are also di- 
agonal. Moreover, we do not take into account the inter- 
band impurity scattering. It is known that the mixing ef- 
fect between different bands disappears both in the weak 
and the strong-coupling limit of the interband impurity 
scattering, hence the multiband situation is restorec&S. 
Therefore, the coupling between bands only appears in 
the gap equation as 

{V- l ) u , *e>R = *TNoiY,(<Ph(k)fai(k; wj) , 

V n 

with A^jj(fe) = A£ij(y3£ 7 (fc), which corresponds to 
eq. <|3*TJl . Here we have used 

= - (V)«, $>MfeM 7 (fe'), (72) 



where the coupling constant V is a matrix in the band 
space. 

Now we express a set of formulas in terms of dimen- 
sionless parameters. The free energy scaled by VoT 2 , 
No = J2e Nqi being the total DOS in the normal state, 
is given by 



9. 
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7V T 2 8tt 



E 



\dt\ 2 \nt + 2nt^2 



\d, 



(73) 
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with 



and 



a ee' = 1 )tt' ~ n t A c 8 u >, 



(74) 



where A = N V, n e = N oe /N , (b, h) = (B, H)/H c2 (T = 
0) and {t,de,o n ,i e ,i' e ) = {T,A e ,w n ,I e ,I' e )/T c (H = 0). 
Since all the argument in the previous section holds for 
each band t, we have added the subscript I to the rele- 
vant quantities. Note that the angular average should be 
carried out over the corresponding Fermi surface. Here 
the quantity defined as k = H C 2/T c ^/Nq is the GL pa- 
rameter, which is of the order of the ratio between the 
penetration depth and the coherence length of the pri- 
mary order parameter (specified by t = 1). The fac- 
tor A c = ln(2e 7 o> c /7rT c ) is the smallest eigenvalue deter- 
mined by 



det a 



(o) 



0. 



(75) 



ft is useful to introduce a dimensionless critical field 
b C 2 = \e\(v/T c ) 2 H C 2/2, which will be determined such 
that the gap vanishes at b = 1 and t — 0. We summa- 
rize appropriate formulas for (a) non s-wave pairings, (b) 
s-wave singlet, and (c) the clean limit. 



(a) non s-wave pairings, 



3i = 



f + 



Vi 



v e ±(k) 



bb, 



(n) , . Vl 

Ve > = ! + -■■ 



V'2 



(ge 



-1/2 



and 



2 9i r- 
U = T~ ; V7T 



o n cos 2 6t + (ge) 2 sin 2 5t ' 
wi \ \d e (ky 2 



(76) 



(77) 



-W(iU£ n 



ve±(k) J \fbb c2 

+ 2o n (4 n) -l)(l- gt ), (78) 



sin Of 



(b) s-wave singlet 



9e = 



i \vt±(k) I bb c2 



(n) , . vt i 

Vt = 1 + — {9e 

On 



(a) 

m 



ge\ T— 777 I W(iue n ) 



sbbrt \ \ve±{k) 



(79) 

-1/2 

1 

(80) 
(81) 

-1 

? 

(82) 



it = 



\ve±{k)J 



\def 



/bb, 



c2 



:W(iue n ) 



geVe 



(a) 



1 + .% 



Vt 



(a) 



+ 2 „(^ n) -i)(i- gt ), (83) 



If = 



o 2 n ( V ^-i) 2 +d 2 (i-i/vn 



(«)\2 



vt 

(c) clean limit 



(84) 



ge = 



i> \ve±(k) I bb C 2 



-1/2 



and 



(85) 



2g e r vi \ \de{k)\ 2 ^ rr 

l-K I — I ^^WilUtn) (86) 



l+ge \v£±(k) I Vbb, 



c2 



where we have defined the impurity scattering rate as 
v t = l/2r m T c and u ln = [vi/v£±(k)]d n /^bb c2 . 

In order to discuss properties under the rotating mag- 
netic field, we move to the new coordinate (x',y',z') in 
which the applied magnetic field is parallel to z' axis. 
Representing the rotation matrix R for the original mag- 
netic field H = H(sm6hCOS(J)h,sm6hSm<f)i l ,cos6f l ) as, 



cos 6h cos fa cos Oh sin fa — sin 6h 
R = ! -smfa cos fa | , (87) 

sin 9h cos fa sin 6h sin fa cos 6h 

we obtain the velocity in the new coordinate as v' e (k) = 
Rvt(k). In the case of H || a; for example, we have 
v 2 , = v 2 and v 2 , = v 2 with 9h = n/2 and fa = as they 
should. Then, we can express the necessary quantities in 
terms of v' e (k) as 

Xt = {vl>(k))/{v%,(k)), (88) 
vt±(k) = yfx7 1/2 vUk) + X 1 e /2 v? v ,(k). (89) 

For details of vt(k) we can refer to band structure calcu- 
lations or simple tight-binding models. 

Using the argument similar to the single-band case, we 
obtain the free energy in the clean limit, 



SN 



y T1 - -h) 2 + J2 (<$ + <*t(t, b)8ee) d\ ■ d t > 

+ |E^( < ' 5 )(l^l 2 ) 2 ' ( 9 °) 
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where 

ae(t, b) = n e 

and 



x W(iu£„) 



(91) 



Pt{t,b) 



-niVK 



•3/2 



n=0 \ V^-L( fc )/ 

x W(iu in )WW(iu in )\ip i7 (k)\ 4 



(92) 



Note that only the quadratic term contains the coupling 
between different bands, which corresponds to the in- 
ternal Josephson coupling. We obtain the temperature 
dependence of the upper critical field, b(t) by solving the 
equation, 



dot (a<$ +ai(t,b)6 u >) = 0. 
In particular, at T — we have 
a t (p,b) = 



ng ( In 




and 



0t(O,b) = n t - 



Vi 



'466 c2 \ \v t± (k) 
Similarly, we obtain the entropy as 



l¥*r(*)| S 



\<P*yW\ 4 



S(t,h) 
N T C 



N{x\t,h) 



-4 / dx 

la N Q 



/(?M? 



with the total DOS, 
N(x;t,h) 



s-waves. For simplicity, we consider strongly type-II 
superconductors (« 3> 1) in the clean limit with two- 
dimensional cylindrical Fermi surfaces. The magnetic 
field is applied along z axis. We also use the same Fermi 
velocity for two bands. Thus we have B = H , x = 1; 
vi/vt±(k) = 1, (• • • ) = f 2n t£ ■■ ■ , and only the gap(s) 
are to be determined. 

Let us first show results for the single-band isotropic 
s-wave, ip(k) = 1. Figure ^ shows the averaged ther- 
modynamic quantities in the mixed state, (a) the field 
dependence of the gap magnitude A scaled by A = 
A(h = t = 0) at different temperatures, (b) the field 
dependence of the DOS at t = 0, (c) the field depen- 
dence of the zero-energy (ZE) DOS at t = 0.1, and (d) 
the temperature dependence of the specific heat normal- 
ized to that of the normal state, Cn ■ In Fig. OJa) the 
gap at the low temperature first increases as the field 
increases. This is an drawback of the present approxi- 
mation, which is not valid at very low fields at low tem- 
perature. The thin line represents an empirical formula, 
A(t,h) = A(t, 0)\/l ~ h, which well describes the field 
dependence of the gap for t > 0.5. In Fig.^b) the peak 
position of the DOS moves upward in energy with the 
external field. In Fig. [TJc) the open circles are taken 
from the results of the full numerical solution of the qua- 
siclassical equations done by Miranovic et alM^l. It 
indicates that the BPT approximation works very well 
for h > 0.6 at low temperatures. For lower fields the 
BPT approximation overestimates contribution from vor- 
tex cores since the vortex core size becomes unphysically 
large in the Abrikosov lattice model, Q37f) (see numerical 
results^ 3 .). 

Next we move to the case of d x 2_ y 2-w&ve, tp(k) = 
(95) VScos^t/)). Figure |U shows results for d x 2_ y 2-wa,ve in 
the plot similar to the previous case. The same empirical 
formula works fine for t > 0.5. A tendency of the peak 
shift in the DOS is the similar to s-wave, but an amount 
of the shift is smaller than that of s-wave. Although the 
BPT approximation overestimates the vortex core con- 
tribution in comparison with full numerical results 54 ' 55 , 
it shows better agreement than the case of s-wave. 

Figure |3 shows results for two-band s-waves. We use 
the same density of states for both bands, m — n 2 = 0.5 
for simplicity. Here the coupling matrix is given by 



(93) 



(94) 



(96) 



iVn 



^ riiRe (gt(k; io n — > x + id) 



(97) 



and the specific heat is given by C(t, h) = t(dS/dt). For 
reference, we have Stf(t,h) — C^(t,h) = 2ir 2 N T/3 in 
the normal state. 



V. EXAMPLES 

Now we demonstrate the application of our theory to 
typical pairings, i.e. (a) the single-band isotropic s-wave, 
(b) the single-band dj.2_, y 2-wave and (c) the two-band 



Ai A 
A A 2 



(98) 



with Ai = 0.25, A 2 = O.O8A1 and A = 0.26Ai. In Fig.^c) 
the ZEDOS's are compared with numerical calculation of 
the Bogoliubov de Gennes framework at T = Gli. The 
field dependence in the passive band agrees very well with 
the full numerical result since the vortex core size is suf- 
ficiently large due to the smallness of the gap with large 
coherence length. The discrepancy mainly comes from 
the primary band but its contribution accounts for ri\ 
of the total ZEDOS. Therefore, the discrepancy in total 
becomes less remarkable. It should be emphasized that 
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behaviors of two-band systems would be changed sensi- 
tively by a slight change of material parameters and a 
combination of two gaps. Therefore in discussing exper- 
imental data it is necessary to take into account precise 
material parameters in accordance with band-structure 
calculations^^. 

Finally, we demonstrate the magnetization curves for 
several k's in the case of s-wave in Fig. 0] k c is the 
critical value of type I-II superconductivity, and n c = 
27r 3 / 2 e- 7 = 6.25278 in our definition of n. The thin lines 
denote results obtained by means of the GL theory^. 
For larger k, the BPT approximation shows better agree- 
ment with the GL results. In the BPT approximation, 
the transition at the lower critical field, H c \ becomes first 
order for k/k c < 1.9, i.e., the free energy has two minima 
in (B,A) space above H c \- In this case, the hysteresis 
loop is represented by the arrow. For small k nonlocal ef- 
fect becomes important. Although the BPT approxima- 
tion properly takes into account nonlocal effect beyond 
the GL theory, the treatment of averaged magnetic field 
itself becomes poor. To conclude the validity of the BPT 
approximation for small k, we should await results of full 
numerical calculation. 



VI. SUMMARY 

We have presented a simple calculational scheme of 
thermodynamic quantities for singlet and unitary triplet 
states under magnetic fields. A combination of the ap- 
proximate analytic solution with a free energy functional 
in the quasiclassical theory provides a wide use formal- 



ism including the impurity scattering and the multiband 
superconductivity. We have discussed the simple formula 
for the upper critical fields in terms of the microscopic 
free energy under magnetic fields. The theory requires 
a little numerical computation as easy as the usual BCS 
theory without magnetic fields. We have demonstrated 
the application to s-wave, d x i _ y i -wave and two-band s- 
wave cases. The comparisons with reliable numerical cal- 
culations conclude that the BPT approximation works 
better for gap with line of nodes than full gap. It also 
works well for multiband superconductivity because of 
its smallness of the gap magnitude in passive band. The 
flexibility of the theory allows us to plug in detail struc- 
ture of the Fermi surface with help of band calculations 
in an arbitrary direction of the magnetic field. All the 
feature of the present theory is appropriate to analyze 
experimental data semi-quantitatively, especially in the 
oriented magnetic field. The method is also useful to dis- 
cuss linear response and transport coefficients, which will 
be discussed in a future publication. 
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FIG. 1: The averaged thermodynamic quantities for the s-wave: (a) the gap magnitude scaled by Ao = A(T = H = 0). The 
thin line represents the empirical formula A(T,H) = A(T, 0)^1 — H/H C 2, (b) the density of states at T = 0, (c) the field 
dependence of the zero-energy DOS compared with full numerical results 53,5 *, (d) the specific heat. 




FIG. 2: The averaged thermodynamic quantities for the d I 2_ y 2-wave: (a) the gap magnitude. Thin line is the same empirical 
formula in Fig. Q (b) the DOS at T = 0, (c) the ZEDOS and full numerical results 5 ^, (d) the specific heat. 
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FIG. 3: The averaged thermodynamic quantities for the two-band s-wave: (a) the gap magnitude scaled by the primary gap 
A = Ai(T = H = 0), (b) the DOS at T = 0, (c) the ZEDOS and numerical results±i, (d) the total specific heat. The coupling 
matrix is given by A = ((Ai, A), (A, A2)) with Ai = 0.25, A2 = O.O8A1 and A = 0.26Ai. 
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FIG. 4: The magnetization curves for several GL parameters. 
The thin lines are taken from the GL result^. For k/k c < 1.9 
the lower transition becomes first order in the BPT approxi- 
mation. The arrow indicates the hysteresis loop. 
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